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Kibble-Zurek mechanism is a theory of defect formation 
in a non-equilibrium continuous phase transition. So far the 
theory has been successfully tested by numerical simulations 
and condensed matter experiments in a number of systems 
with small thermal fluctuations. This paper reports first nu- 
merical test of the mechanism in a system with large thermal 
fluctuations and strongly non-mean-field behavior: the two 
dimensional Ising model. The theory predicts correctly the 
initial density of defects that survive a quench from the disor- 
dered phase. However, before the system leaves the Ginzburg 
regime of large fluctuations most of these defects are annihi- 
lated and the final density is determined by the dynamics of 
the annihilation process only. 
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Introduction. — In a system with a continuous phase 
transition an adiabatic change of a parameter of the sys- 
tem, like e.g. temperature, pressure or a coupling con- 
stant in a Hamiltonian, can drive the system from a 
disordered phase to an ordered one. A classic example 
is the paramagnet-ferromagnet transition in the two di- 
mensional (2D) Ising model. Thermodynamics of con- 
tinuous phase transitions has been intensively explored 
over many years. Two mayor achievements: the solution 
of the Ising model by Onsager and the renormalization 
group of Wilson were rewarded with a Nobel Prize in 
physics. The RG formalism revealed deep connections 
between statistical mechanics and quantum field theory. 

A candidate theory of non-equilibrium phase transi- 
tions is the Kibble-Zurek mechanism (KZM) [1,2]. Kib- 
ble pointed out [1] that in a non-equilibrium transition 
there is no time to fully develop the long range order 
characteristic for the ordered phase. As a result, the fi- 
nal state of the system is a mosaic of finite size ordered 
domains with different orientations of the order param- 
eter in every domain. In a topologically non-trivial case 
this disorder takes the form of a finite density of topo- 
logical defects. This qualitative idea was quantified more 
by Zurek [2] . Zurek mechanism is a combination of two 
basic facts: (1) a divergent correlation length 



(i) 



where e is a dimensionless distance from the critical point, 
v is a critical exponent, and £o is a microscopic length 
scale, and (2) the critical slowing down or divergent re- 
laxation time 



Here To is a microscopic time scale. Because of the diver- 
gent relaxation time any finite rate transition is a non- 
adiabatic phase transition: sufficiently close to the crit- 
ical point (where e = 0) the system is too slow to react 
to the changing external parameter e(t). Close to e = 
we can linearize 



e(t) 



t 
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The relaxation time (2) equals the transition rate |e/^ 
at ez ~ {tq/tq)~ when the correlation length (1) is 



(4) 
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This Zurek length is the average size of the ordered do- 
mains after the phase transition and it determines the ini- 
tial density of topological defects frozen into the ordered 
phase after a non-adiabatic continuous phase transition. 

The original motivation for Kibble and Zurek were 
symmetry breaking phase transitions in cosmology. The 
random topological defects arising in such transitions 
might provide initial seeds for structure formation in the 
early Universe [3]. However, the universality of phase 
transitions makes these ideas also relevant for a wide va- 
riety of condensed matter systems where they can be ver- 
ified by experiment. 

The KZM prediction (4) is supported by a number of 
numerical simulations [4]. However, as a result of finite 
numerical resources these numerical data arc limited to 
fast quenches (small tq) with a large ez in the regime of 
small fluctuations where one can use the mean field (MF) 
value of the critical exponent %p = i . KZM is also sup- 
ported by experiments in systems with small fluctuations 
like superfluid helium 3 [5], low T c superconductors [6], 
and fast quenches in high T c superconductors [7]. In 
contrast, experiments in systems with large fluctuations 
like helium 4 [8] or slow quenches in high T c supercon- 
ductors [9] are inconclusive. Rivers suggested [10] that 
vortices created in the helium 4 experiment [8] disappear 
in a faster than expected annihilation. Due to technical 
difficulties the analytic calculations in Ref. [10] eventu- 
ally resort to a linearization equivalent to the mean-field 
theory. It is suggested there that beyond this linearized 
theory close to the critical point the annihilation rate is 
divergent. However, simulations in Ref. [11] show that 
this effect may be not as dramatic as anticipated in Ref. 
[10]. These authors suggest that because of the critical 
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slowing down the annihilation rate close to e = may 
in fact vanish. Due to limited numerical resources the 
numerical evidence in Ref. [11] is rather indirect. To 
summarize, the problem of KZM in the Ginzburg regime 
of large fluctuations has been recognized [10] but is far 
from being settled. 

At the moment we do not have any condensed matter 
or numerical experiment supporting KZM for large fluc- 
tuations and at the same time this is the regime where 
KZM in principle can be questioned on general grounds. 
The argument leading to Eq.(4) implicitly assumes that 
close to the critical point the divergent correlation length 
£ in Eq.(l) is the only relevant length scale. However, as 
is well known [12] but not quite generally appreciated, if £ 
were the only length scale, then, on dimensional grounds, 
all critical exponents would take their mean field values. 
As they do not (for example, in the 2D Ising model v — 1 
instead of the mean field ^mf = 5), then both £ and 
the microscopic £0 must be relevant. With two relevant 
length scales the dimensional argument alone is not suf- 
ficient to determine the initial density of defects. 

In this paper I report first numerical test of KZM for 
large fluctuations. As the critical regime is numerically 
demanding (large £ means large lattice and large r means 
long time) I chose the simplest possible model - the cel- 
ebrated 2D Ising model. This simple model has v — 1 
clearly different from the mean field i^mf = 1/2, and it 
has no regime where the MF theory might be at least re- 
motely accurate. It is a perfect testing ground for KZM. 

Ising model with Glauber dynamics. — Hamilto- 
nian of the ferromagnetic Ising model in 2D is 

H = ]T S, .S ; . (5) 

{id) 

Spins Si e {— 1, +1} sit on a 2D N x N lattice with peri- 
odic boundary conditions, means a pair of nearest 
neighbor sites. The microscopic lengthscale £ = 1 is the 
lattice spacing. In all the following numerical simulations 
a 1024 x 1024 lattice was used. 

To study non-equilibrium phase transitions the Ising 
model has to be supplemented with dynamics. The stan- 
dard choice is Glauber dynamics also known as Monte- 
Carlo with a heat bath [13]. In the Glauber algorithm 
every time step consists of the following sub-steps: 

• choose a random spin Si from the lattice, 

• calculate its local field hi — — J2j n n % &j , 

• calculate a probability P = exp(/3/ij), 

• choose a random number r € [0, 1], 

• if r > P then set Si = +1, else set S 1 , = — 1. 

Here (3 is an inverse temperature of the heat bath. This 
algorithm relaxes the state of the Ising model towards 



thermal equilibrium at a given (3 [13]. On average it 
takes A 2 steps to update the state of N 2 spins on the 
lattice. These A 2 steps define the microscopic time scale 
To which I set equal to 1. 

The Ising model with Glauber dynamics belongs to the 
same universality class as the (j) 4 model with noise r\ 

to J^ = £ o 2 V 2 0-A(0 2 - 1)^ + 7? (6) 

so often employed in the numerical simulations of KZM 
[4]. Here the continuum real field <j> is a coarse grained 
lattice spin Si. The Ising model is an efficient "molecular 
dynamics" version of the </> 4 field theory (6). 

Relaxation time. — In order to estimate the expo- 
nent y in Eq.(2) the relaxation time r was measured for 
several values of (3 < (3 C . For each (3 the Ising model was 
initially prepared in a fully polarized state with all Si — 
1, and then its average magnetization M = ^2iSi/N 2 
relaxed towards the equilibrium at M = 0, see the insert 
in Fig.l. Each magnetization decay was fitted with an 
exponent M = cxp(— t/r). The best fits of r are shown 
in the double logarithmic Fig.l as a function of (3 C — [3. 
The slope of the linear fit in Fig.l gives an estimate of 
y = 2.09 ±0.02. 
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FIG. 1. y = log 10 r as a function of x = log 10 (/3 c — /?). The 
rs are the best fits to the exponential decays of magnetization 
shown in the insert. The solid line is the best linear fit with 
a slope of y = 2.09 ± 0.02. 

Quenches. — Phase transitions were simulated with 
a linear ramp of the inverse temperature 

(3{t) = 1.5 — . (7) 

The initial state at t = was a state with random mu- 
tually uncorrelated spins - the state of equilibrium at 
(3 = 0. Fig. 2 shows density of domain walls separating 
positive Si from negative Si as a function of (3 for a num- 
ber of different transition times tq. The critical point 
is (3 C = 0.4407. For large tq the density plots approach 
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the equilibrium density n eq ((3). Note that the equilib- 
rium density n eq ((3) of domain walls remains nonzero 
even for (3 > [3 C . This is the critical Ginzburg regime 
of large fluctuations. A non-equilibrium transition with 
a finite tq results in an additional non-equilibrium den- 
sity dn(/3) = n(/3)- n eq {(3) > 0. KZM predicts that 



dn 



KZM 



(/?) 



HI 



-0.324±0.003 



(8) 



Before this prediction is compared with the numerical 
data in Fig. 2, let me digress on annihilation of domain 
walls. 
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FIG. 2. Total density n(/3) of domain walls separating pos- 
itive and negative Si as a function of /3 for several values of 
the quench time tq = 2, 4, 8, ... , 65536 (from top to bottom). 
For the initial state of random spins the density is normalized 
to n = 1. For large tq the plots tend to the equilibrium den- 
sity of defects n eq (/3) which is finite even for j3 > f3 c — 0.4407 
in the Ginzburg regime of large fluctuations. 

Defect annihilation. — First example is annihilation 
of defects from an initially totally random spin config- 
uration. The initial dn(t = 0) decays in time. Fig. 3 
shows the equilibrating n{t) for several values of (3 > (3 C . 
Each decay is fitted with a solid line that follows the 
power law dn(t) — (Ta/t) 1 / 2 with an exponent of 1/2 
known from the theory of phase ordering kinetics [14]. 
The best fits are r a = 0.86 ± 0.05, 0.93 ± 0.05, 0.64 ± 0.05 
for (3 = 0.47, 0.60, 1.0 respectively. They are more or less 
constant in the considered range of temperatures: as the 
critical point is approached the time scale for annihila- 
tion T a neither diverges (as suggested in Ref. [10]) nor 
vanishes (as suggested in Ref. [11]), but remains finite 
and close to the microscopic tq = 1, 



TO 



(9) 



The quench time tq determines the time available for 
defect annihilation. At late times after the transition, 
when most of the initial KZM domain walls are already 
annihilated, we expect the scaling 



dn{f3) 



(10) 



It also follows from a phcnomcnological equation: 
-7}dn 3 (t). Its solution is 



T £dn(t) 



dn(t) = 



dn{0) 



(11) 




Note that at late times dn(t) is forgetting its initial value 
dn(0) = driKZM- This solution is an ilustration of the 
exact result (10) from Ref. [14]. 

Second example is annihilation of domain walls from 
an initial state of equilibrium at (3 > (3 C . The ini- 
tial state was prepared by starting from fully polar- 
ized spins, all Sj = 1, and then heating them up at 
(3 = 0.45 for a time of 10 5 sufficient to reach thermal 
equilibrium with n eg (0.45) = 0.20. Then at t = (3 
was suddenly increased (the heat bath was cooled) to 
(3 = 0.55. Fig. 3 shows n(t) decaying towards the new 
equilibrium at n eq (0.55) — 0.075. This decay is much 
faster than for random initial spins because the equilib- 
rium domain walls in the Ginzburg regime at (3 > (3 C 
are just boundaries of bubbles of the minority spin-down 
phase in the spin-up polarized fcrromagnet. The bubbles 
together with their walls decay soon after the tempera- 
ture is turned down. 
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FIG. 3. Density of defects n(t) starting from an initial 
state with random spins and decaying towards n eq (/3) for 
/3 = 0.47,0.6, 1.0. The "0.45 ~* 0.55" marks the plot of n(t) 
starting from the state of equilibrium at (3 = 0.45 > (3 C in the 
Ginzburg regime and decaying quickly towards a new equilib- 
rium at (3 = 0.55. 

KZM versus annihilation. — Figure 4 is a double 
logarithmic plot of the non-equilibrium density dn{(3) in 
Fig. 2 as a function of tq for a number of (3s. The slope 
at the critical (3 C = 0.4407 is -0.315 ± 0.007. This slope 
is consistent with the KZM slope (8) of —0.324 and very 
different from a mean-field KZM slope of —0.65 for cmf = 
1/2. The initial non-equilibrium density of domain walls 
is determined by KZM. 

In contrast, similar slopes for (3 — 1.0 and 1.5 are 
—0.45 ± 0.01 and -0.48 ± 0.01 respectively, and they are 
consistent with the phase ordering kinetics exponent of 
— 1/2 in Eq.(10). Apparently at later times the system 
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forgets the initial density driKZM and dn((3) is determined 
solely by the dynamics of defect annihilation. 

Indeed, the circles in Fig. 4 show dn((3 — 1.5) for a 
simulation where (3{t) is ramped up like in Eq.(7), but 
starting from the initial /?o = 0.6 > /3 C . The spins were 
random at the initial /3 - The circles sit on the solid line 
which is a fit to dn(/3 — 1.5) obtained from a full quench 
like in Eq.(7). The annihilation dominated dn{j3) at later 
times is not sensitive to the details of the KZM of defect 
formation, compare Eqs.(10,ll). 

However, the defects that survive annihilation at later 
times are KZM defects quenched in from the disordered 
phase. As we have already seen, compare Fig. 3, that an- 
nihilation of the Ginzburg domain walls is much faster 
than annihilation of defects from the initially random 
spin state. The latter state contains large domain walls 
while in the former domain walls are boundaries of bub- 
bles of a minority spin phase. The points in Fig. (4) con- 
nected by a dashed line show dn([3 = 1.5) after a quench 
starting from the equilibrium state at /3q = 0.47 in the 
Ginzburg regime. These densities are orders of magni- 
tude lower than densities from the full quenches starting 
at j3 = 0: Ginzburg defects do not survive annihilation. 



fects is gradually annihilated and when the system leaves 
the Ginzburg regime the density of defects is no longer 
sensitive to the details of the KZM, but it is determined 
by the dynamics of the annihilation process only. In 
particular, the dependence of the density on the transi- 
tion rate is determined by an exponent that comes from 
the theory of phase ordering kinetics and not from the 
KZM. The only way to see the KZM scaling (8) directly 
is to measure the amount of disorder close to the critical 
point where the non-equilibrium KZM density is largely 
obscured by the prevailing equilibrium thermal fluctua- 
tions. However, the defects that survive the annihilation 
are the KZM defects quenched in from the high temper- 
ature phase, the defects quenched in from the Ginzburg 
regime decay much faster. The surviving defects are a 
clear, though indirect, signature of the KZM. 
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FIG. 4. y — log 10 dn(/3) as a function of x = log 10 tq 
for f3 = 0.4407, 1.0, 1.5 from top to bottom. 
Solid lines are the best linear fits with slopes of 
-0.315 ± 0.007, -0.45 ± 0.01, -0.48 ± 0.01 respectively. Cir- 
cles show dn((5 = 1.5) in a quench starting from f3o = 0.6 
and random initial spins. The points connected by a dashed 
line show densities dn(/3 — 1.5) in a quench starting from 
/3o = 0.47 in the Ginzburg regime and spins initially in ther- 
mal equilibrium. 

Conclusion. — I presented first numerical test of the 
Kibble-Zurek mechanism (KZM) in the Ginzburg regime 
of large thermal fluctuations. In this regime both the 
Zurek length £z and the microscopic length £o are rele- 
vant length scales that determine the density of defects. 
However, the density of non-equilibrium defects frozen 
into the ordered phase by a quench from the disordered 
phase is determined by £z only. This initial density of de- 
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